• Session 1: Motivation, why and how to think about data, and getting started with R
• Session 2: Making basic plots, grammar of graphics, good practices
• Session 3: Advanced graphics, layering, using maps

(If you re-started RStudio, be sure to re-open your project too.)

Read the OECD PISA data

student2012.sub <- readRDS("../data/student_sub.rds")
dim(student2012.sub)
#> [1] 271323     50
student2012.sub$ST04Q01 <- factor(student2012.sub$ST04Q01, 
  levels=c(1,2), labels=c("Female", "Male"))

Calculate the statistics: Mean difference

student2012.stats <- student2012.sub %>% 
  group_by(CNT) %>%
  summarise(wmathgap=weighted.mean(PV1MATH[ST04Q01=="Male"], 
                  w=SENWGT_STU[ST04Q01=="Male"], na.rm=T)-
               weighted.mean(PV1MATH[ST04Q01=="Female"],
                  w=SENWGT_STU[ST04Q01=="Female"], na.rm=T))
kable(head(student2012.stats))

• Means are point estimates
• Map mean difference to point geom, using position along a line

ggplot(data=student2012.stats) + 
  geom_point(aes(x=CNT, y=wmathgap), size=3) + 
 coord_flip() + theme_bw()

• Arrange: by country
• Order by magnitude of difference

student2012.stats$CNT <- factor(student2012.stats$CNT, 
    levels=student2012.stats$CNT[order(student2012.stats$wmathgap)])
ggplot(data=student2012.stats) + 
  geom_point(aes(x=CNT, y=wmathgap), size=3) + 
 coord_flip() + theme_bw()

Create bootstrap confidence intervals for each mean difference

cifn <- function(d, i) {
  x <- d[i,]
  ci <- weighted.mean(x$PV1MATH[x$ST04Q01=="Male"], 
          w=x$SENWGT_STU[x$ST04Q01=="Male"], na.rm=T)-
        weighted.mean(x$PV1MATH[x$ST04Q01=="Female"],
          w=x$SENWGT_STU[x$ST04Q01=="Female"], na.rm=T)
  ci
}

bootfn <- function(d) {
  r <- boot(d, statistic=cifn, R=100)
  l <- sort(r$t)[5]
  u <- sort(r$t)[95]
  ci <- c(l, u)
  return(ci)
}

Apply ci functions to data

student2012.sub.summary.gap.boot <- student2012.sub %>% 
  split(.$CNT) %>% purrr::map(bootfn) %>% data.frame() %>%
  gather(CNT, value)
student2012.sub.summary.gap.boot$ci <- 
  rep(c("ml","mu"), 
      length(unique(student2012.sub.summary.gap.boot$CNT)))
student2012.sub.summary.gap.boot.wide <- 
  student2012.sub.summary.gap.boot %>% 
  spread(ci, value)
student2012.sub.summary.gap <- merge(student2012.stats,
  student2012.sub.summary.gap.boot.wide)
kable(head(student2012.sub.summary.gap))

ggplot(data=student2012.sub.summary.gap) + 
  geom_point(aes(x=CNT, y=wmathgap), size=3) + 
  geom_segment(aes(x=CNT, xend=CNT, y=ml, yend=mu)) + 
  coord_flip() + theme_bw() 

Match three digit codes to country names, more recognizable labels

student2012.sub.summary.gap$name <- NA
for (i in 1:length(student2012.sub.summary.gap$name))  
  student2012.sub.summary.gap$name[i] <-
  isoToName(as.character(student2012.sub.summary.gap$CNT[i]))
# QCN is Shanghai, not whole of China - 
# Don't know what country TAP is
student2012.sub.summary.gap$name[
  student2012.sub.summary.gap$CNT == "QCN"] <- 
  isoToName("CHN")
student2012.sub.summary.gap$name[
  student2012.sub.summary.gap$CNT == "TAP"] <- 
  "TAP"

Create categorical gap variable to draw attention to significant difference

student2012.sub.summary.gap$wmathgap_cat <- "same"
student2012.sub.summary.gap$wmathgap_cat[
  student2012.sub.summary.gap$ml > 0] <- "boys"
student2012.sub.summary.gap$wmathgap_cat[
  student2012.sub.summary.gap$mu < 0] <- "girls"
kable(head(student2012.sub.summary.gap))

student2012.sub.summary.gap$name <- factor(student2012.sub.summary.gap$name, 
    levels=student2012.sub.summary.gap$name[
      order(student2012.sub.summary.gap$wmathgap)])
kable(head(student2012.sub.summary.gap))

ggplot(data=student2012.sub.summary.gap) + 
  geom_hline(yintercept=0, colour="grey80") + 
  geom_point(aes(x=name, y=wmathgap, color=wmathgap_cat), 
             size=3) + 
  geom_segment(aes(x=name, xend=name, y=ml, yend=mu, 
                   color=wmathgap_cat)) + 
  coord_flip() + theme_bw() 

• Labels
• Axis limits
• Grid lines
• Color

ggplot(data=student2012.sub.summary.gap) + 
  geom_hline(yintercept=0, colour="grey80") + 
  geom_point(aes(x=name, y=wmathgap, color=wmathgap_cat), size=3) + 
  geom_segment(aes(x=name, xend=name, y=ml, yend=mu, 
     color=wmathgap_cat)) + xlab("") +  
  scale_colour_manual("", values=c("boys"="skyblue", 
    "girls"="pink", "same"="lightgreen")) +
  scale_y_continuous("Girls <----------> Boys", 
    breaks=seq(-30, 30, 10), limits=c(-35, 35), 
    labels=c(seq(30, 0, -10), seq(10, 30, 10))) + 
  coord_flip() + theme_bw() + 
  theme(axis.text.x = element_text(size=5), 
        axis.text.y = element_text(size=5), 
        axis.title = element_text(size=7), 
        legend.text = element_text(size=5),
        legend.title = element_text(size=5))

ggplotly()

Map data is essentially a set of points, and line segments. You can get maps from various sources, and wrangle the files/data into an R object. This can be merged with data to provide spatial context to problems.

world <- getMap(resolution = "low")
extractPolys <- function(p) {
  polys <- NULL
  for (i in 1:length(p)) {
    for (j in 1:length(p[[i]]@Polygons)) {
      x <- p[[i]]@Polygons[[j]]@coords
      polys$lon <- c(polys$lon, x[,1])
      polys$lat <- c(polys$lat, x[,2])
      polys$ID <- c(polys$ID, rep(p[[i]]@ID, nrow(x)))
      polys$region <- c(polys$region, 
        rep(paste(p[[i]]@ID, j, sep="_"), nrow(x)))
      polys$order <- c(polys$order, 1:nrow(x))
    }
  }
  return(data.frame(polys))
}
polys <- extractPolys(world@polygons)

Here is what is looks like:

kable(head(polys))

ggplot(data=filter(polys, region=="China_1"), 
       aes(x=lon, y=lat)) + 
  geom_point()

ggplot(data=filter(polys, region=="China_1"), 
       aes(x=lon, y=lat, order=order)) + 
  geom_path()

ggplot(data=filter(polys, region=="China_1" | region=="Australia_1"), 
       aes(x=lon, y=lat, group=region, order=order)) + 
  geom_path()

ggplot(data=polys) + 
  geom_path(aes(x=lon, y=lat, group=region, order=order))

Join education data with map polygons

polys <- polys %>% rename(name = ID)
student2012.sub.map <- left_join(
  student2012.sub.summary.gap, polys)
student2012.sub.map <- student2012.sub.map %>% 
  arrange(region, order)

Make it look like a map, by tweaking the plot appearance

theme_map <- theme_bw()
theme_map$line <- element_blank()
theme_map$strip.text <- element_blank()
theme_map$axis.text <- element_blank()
theme_map$plot.title <- element_blank()
theme_map$axis.title <- element_blank()
theme_map$panel.border <- element_rect(
  colour = "grey90", size=1, fill=NA)

ggplot(data=polys) + 
  geom_path(aes(x=lon, y=lat, group=region, order=order), 
            colour=I("grey90"), size=0.1) + 
  geom_polygon(data=student2012.sub.map, aes(x=lon, y=lat, 
            group=region, order=order,  
            fill=wmathgap_cat)) +
  scale_fill_manual("Diff>5", values=c("boys"="skyblue", 
                                    "girls"="pink", 
                                    "same"="lightgreen")) + 
  scale_x_continuous(expand=c(0,0)) + 
  scale_y_continuous(expand=c(0,0)) +
  coord_equal() + theme_map 

ggplotly()

Image maps can also be pulled using the package ggmap to use as the background for data.

gm <- get_googlemap(center = c(lon=144.96, lat=-37.815), zoom=14)
ggmap(gm) 

sensor_loc <- read_csv("../data/Pedestrian_Sensor_Locations.csv")
ggmap(gm) + geom_point(data=sensor_loc, aes(x=Longitude, y=Latitude))

• With your neighbor, take one of the problems that you brainstormed on the first day, and work a question through from start to end, summarise.
• Choices were: economics, gapminder, pedestrian sensor.

melbtemp <- read.fwf("../data/ASN00086282.dly", 
   c(11, 4, 2, 4, rep(c(5, 1, 1, 1), 31)), fill=T)
melbtemp.m <- melbtemp %>%
  select(num_range("V", c(1,2,3,4,seq(5,128,4)))) %>%
  filter(V4 %in% c("PRCP", "TMAX", "TMIN")) %>%
  gather(day, value, V5:V125, na.rm = TRUE) %>%
  spread(V4, value) %>%
  mutate(
    tmin = as.numeric(TMIN) / 10,
    tmax = as.numeric(TMAX) / 10,
    t_range = tmax - tmin,
    prcp = as.numeric(PRCP) / 10
  ) %>%
  rename(stn=V1, year=V2, month=V3)

melbtemp.m$day <- factor(melbtemp.m$day, 
  levels=c("V5","V9","V13","V17","V21","V25","V29",
           "V33","V37","V41","V45","V49","V53","V57",
           "V61","V65","V69","V73","V77","V81","V85",
           "V89","V93","V97","V101","V105","V109",
           "V113","V117","V121","V125"),
  labels=1:31)
melbtemp.m$date <- as.Date(paste(melbtemp.m$day, 
     melbtemp.m$month, melbtemp.m$year, sep="-"),
     "%d-%m-%Y")

What are the temperature patterns in Melbourne over time? Daily resolution: global and seasonal patterns

kable(head(melbtemp.m))

ggplot(data=melbtemp.m, aes(x=date, y=tmax)) + geom_point(size=0.1)

melbtemp.m <- melbtemp %>%
  select(num_range("V", c(1,2,3,4,seq(5,128,4)))) %>%
  filter(V4 %in% c("PRCP", "TMAX", "TMIN")) %>%
  gather(day, value, V5:V125, na.rm = TRUE) %>%
  mutate(value=ifelse(value<(-9000), NA, value)) %>%
  spread(V4, value) %>%
  mutate(
    tmin = as.numeric(TMIN) / 10,
    tmax = as.numeric(TMAX) / 10,
    t_range = tmax - tmin,
    prcp = as.numeric(PRCP) / 10
  ) %>%
  rename(stn=V1, year=V2, month=V3)

ggplot(data=melbtemp.m, aes(x=date, y=tmax)) + geom_point(size=0.1)

Aspect ratio emphasizes global trend

ggplot(data=melbtemp.m, aes(x=date, y=tmax)) + 
  geom_point(size=0.5, alpha=0.4) +
  theme_bw() + theme(aspect.ratio=2)

Aspect ratio emphasizes seasonality

ggplot(data=melbtemp.m, aes(x=date, y=tmax)) + 
  geom_point(size=0.1) +
  theme(aspect.ratio=0.1)

ggplot(data=melbtemp.m, aes(x=factor(month), y=tmax)) + 
  geom_boxplot() +
  theme(aspect.ratio=0.3)

Cooler in June-July. Variability is higher in the summer.

Long time series make it difficult to examine trend. Aggregate on year, and use the min and max values.

melbtemp.m.r <- melbtemp.m %>% 
  group_by(year) %>%
  summarise(tmin=min(tmin, na.rm=TRUE), 
            tmax=max(tmax, na.rm=TRUE)) 

ggplot(data=melbtemp.m.r, aes(x=year, ymin=tmin, ymax=tmax)) + 
  geom_ribbon(alpha=0.5) + 
  geom_smooth(aes(x=year, y=tmax), se=F, method="lm", colour="black") +
  geom_smooth(aes(x=year, y=tmin), se=F, method="lm", colour="black") +
  theme(aspect.ratio=1) + xlab("") + ylab("Temperature")

ggplot(data=filter(melbtemp.m, year > 1990, year < 2011), 
       aes(x=factor(month), y=tmax)) + 
  geom_boxplot() + facet_wrap(~year) +
  theme(aspect.ratio=0.3)

Examine global and seasonal patterns of precipitation in Melbourne

Occasionally you would like to organize your plots in special ways. The gridExtra can be used to take individual plots and lay them out together.

p1 <- ggplot(data=subset(CO2.all, flg < 2), 
             aes(x=date, y=co2, colour=stn)) +
  geom_line() + xlab("Year") + ylab("CO2 (ppm)") + 
        facet_wrap(~stn, ncol=1) + 
  theme(axis.text.y=element_text(size = 6), legend.position="none")
p2 <- ggplot(data=subset(CO2.all, flg < 2), 
             aes(date, co2, colour=stn)) +
  geom_line() + xlab("Year") + ylab("CO2 (ppm)") + 
  theme(axis.text.y=element_text(size = 6), legend.position="none")
p3 <- ggplot(data=polys) + 
  geom_path(aes(x=lon, y=lat, group=region, order=order), 
            colour=I("grey60"), size=0.1) + 
  geom_point(data=CO2.all.loc, aes(x=lon, y=lat, group=1), colour="red", 
                      size=2, alpha=0) +
  geom_text(data=CO2.all.loc, aes(x=lon, y=lat, label=stn, group=1), 
            colour="orange", size=5) +
  coord_equal() + theme_map 

grid.arrange(p1, p2, p3, layout_matrix = rbind(c(1,2),c(3,3)))

Notes prepared by Di Cook, building on joint workshops with Carson Sievert, Heike Hofmann, Eric Hare, Hadley Wickham.

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.